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Generalized ABCD propagation for interacting atomic clouds. 
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We present a treatment of the nonlinear matter wave propagation inspired from optical methods, 
which includes interaction effects within the atom optics equivalent of the aberrationless approxi- 
mation. The atom-optical ABCD matrix formalism, considered so far for non-interacting clouds, is 
extended perturbatively beyond the linear regime of propagation. This approach, applied to dis- 
cuss the stability of a matter-wave resonator involving a free-falling sample, agrees very well with 
the predictions of the full nonlinear paraxial wave equation. An alternative optical treatment of 
interaction effects, based on the aberrationless approximation and suitable for cylindrical paraxial 
beams of uniform linear density, is also adapted for matter waves. 

PACS numbers: 03.75.Pp ,03.75.-b, 42.65.Jx, 41.85.Ew, 31.15-Md 
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I. INTRODUCTION 

Light and matter fields are governed by similar equa- 
tions of motion Both photons and atoms interact 
in a symmetrical manner: atom-atom interactions are 
mediated through photons, while photon-photon inter- 
actions are mediated through atoms. Before the advent 
of Bose-Einstein condensation, two groups realized inde- 
pendently that atomic interactions give rise to a cubic 
nonlinearity in the propagation equation analogous to 
that induced by the Kerr effect Following this anal- 

ogy, the field of non-linear atom optics emerged in the 
last decade, leading to the experimental verification with 
matter waves of several well-known nonlinear optical 
phenomena \i\ : the four- wave mixing [9(, the formation 
of solitons ESJS E El and of vortices 0, 0]|, the 
superradiance [l6l ] and the coherent amplification [TtJ • 
The nonlinear propagation of matter waves has been the 
object of extensive experimental [H, [l9[ and theoretical 
work, among which the time-dependent Thomas-Fermi 
approximation [20j . the variational approach pi ], and 
the method of moments [22j . These treatments have 
been used successfully to obtain analytical expressions 
in good agreement with the exact solution of the 3D 
nonlinear Schrodinger equation (NLSE). 

There exists, for cylindrical wave-packets propagat- 
ing in the paraxial regime, a very elegant method to 
handle this equation which has been used in optics 
to treat self- focusing effects [H, H3|- It relies on the 
"aberrationless approximation" , assuming that the 
nonlinearity is sufficiently weak as to preserve the shape 
of a fundamental Gaussian mode, and it involves a 



[i] Many other optical phenomena have also been verified with mat- 
ter waves. A short list includes interferences [5J and diffraction 
phenomena 0, [f|, the temporal Talbot effect [f|, and the in- 
fluence of spatial phase fluctuations on interferometry • New 
effects arise also with rotating condensates |8|. 



generalized complex radius of curvature. This treatment 
is equally relevant for the paraxial propagation of 
cylindrical matter waves, and it is presented in this 
context in Appendix [X] Unfortunately, the assumptions 
required - such as the constant longitudinal velocity 
and the paraxial propagation - limit the scope of 
this approach, which appears as too stringent to de- 
scribe the matter wave propagation in most experiments. 



This motivates the introduction of a different ana- 
lytical method to obtain approximate solutions for the 
NLSE in a more general propagation regime. This is 
the central contribution of this paper, which exposes a 
perturbative matrix analysis especially well-suited to 
discuss the stability of a matter-wave resonator. With 
an Hamiltonian quadratic in position and momentum 
operators, and in the absence of atomic interactions, 
the Schrodinger equation admits a basis of Gaussian 
solutions. Their evolution is easily obtained thr ough 
a time-dependent matrix denoted "ABCD" [J, |25| . 
in analogy with the propagation of optical rays in 
optics [26] . In the "aberrationless approximation" , it 
is possible to extend this treatment to include pertur- 
batively interaction effects and obtain the propagation 
of a fundamental Gaussian mode with a modified 
"ABCD" matrix. As an illustration of this method, 
the stability of a matter-wave resonator is analyzed 
thanks to this "ABCD" matrix, which encapsulates the 
divergence resulting from the mean-field potential. An 
ABCD-matrix approach had already been used in [l8[ 
to characterize the divergence of a weakly outcoupled 
atom laser beam due to interactions with the source 
condensate. The present treatment is sensibly different, 
since it is not restricted to the paraxial regime and 
since it addresses rather self-interaction effects in the 
beam propagation. An "ABCD" matrix, including 
self-focusing effects, is computed in Sec. IIV1 and used 
to model the propagation of an atomic sample in a 
matter-wave resonator. Self-focusing is also discussed 
through an alternative method exposed in Appendix [XJ 
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non-linear evolution: 



Our approach is indeed mainly inspired from pre- 
vious theoretical developments in optics, which aimed 
at treating the wave propagation in a Kerr medium 
through such a matrix formalism [23| . An approach 
of the non-linearity based on the resulting frequency- 
dependent diffraction [27[ successfully explained the 
asymmetric profile of atomic and molecular intra-cavity 
resonances [28|, as well as the dynamics of Gaussian 
modes in ring and two-isotopes lasers [2i|[3(|. Later, a 
second-order polynomial determined by a least-square fit 
of the wave intensity profile was considered to model the 
Kerr effect [3lJ. In this paper, we explore the quantum 
mechanical counterpart of this strategy: mean-field in- 
teractions are modelled thanks to a second-order polyno- 
mial, determined perturbatively from the wave-function, 
and which can be interpreted in optical terms. 



II. LENSING POTENTIAL 

One considers the propagation of a zero-temperature 
condensate in a uniform gravity field and in the mean- 
field approximation. The corresponding Hamiltonian 
reads: 



H = ^-+mgz + g I \(j ) {v 1 t)\ 2 
Am 



(1) 



gi is the coupling constant related to the s-wave scat- 
tering length a and to the number of atoms N by gi — 
ATrNh 2 a/m. Our purpose is to approximate the mean- 
field potential gi\4>(r,t)\ 2 by an operator leading to an 
easily solvable wave equation and as close as possible 
to the interaction potential. A second-order polynomial 
in the position and momentum operators is a suitable 
choice, since it allows to obtain Gaussian solutions to the 
propagation equation. These solutions are approximate, 
but they lead nonetheless to a satisfactory description 
of the propagation of diluted atomic wave-packets and 
of their stability in resonators, which are the issues ad- 



dressed in this paper. We note Ho 
interaction-free Hamiltonian and 

H(r,p,t) =H +P l (r,p,t) 



2 m 



mgz the 



(2) 



E{P{t),\<t>(t))) 



d 3 r | (r|P(f, p, t)\<f>(t)) - gi\4>(T, t)\ 2 j>(r, t)f 

(3) 



The minimization of the distance E ( P(t) , \<fi{t)) ) for 
the lensing potential Pi(t) implies that the function E 
is stationary towards any second-order polynomial coef- 
ficient at the point (Pi(t), \ <f>(t))): 

Vt > t VpEiPit), \<t>(t))) = (4) 

We have noted V p the gradient associated with the coef- 
ficients of a second-order polynomial, and to is the initial 
time from which we compute the evolution of the wave- 
function - we assume that (f>(r,to) is known -. The de- 
termination of the lensing potential associated with self- 
interactions in the beam indeed requires previous knowl- 
edge of the wave-function evolution. This difficulty did 
not arise in other optical treatments of atomic interaction 
effects [H, |H, [33[ , in which the atomic beam propaga- 
tion was mainly affected by interactions with a different 
sample of well-known wave-function. This is typically 
the case for a weakly outcoupled continuous atom laser 
beam, in which the diverging lens effect results from the 
source condensate. We propose to circumvent this self- 
determination problem thanks to a perturbative treat- 
ment. Such approach is legitimate for the diluted matter 
waves involved in usual atom interferometers. The first- 
order lensing polynomial and the corresponding Hamil- 
tonian H^(t) = Hq + P ; (1) (f,p,i) are determined from 
the linear evolution, according to: 

Vt > t V P E (p/ X) (f) , e - l/ ^ o(t - to) |<^ ))) = 

(5) 

Higher-order lensing effects can be computed iteratively. 

(2) 

For instance, the second-order lensing polynomial P ; (t) 
satisfies at any instant t > to'. 



V P P(p/ 2) (i),T 



-i/nf* o dt'iH +p^^,f>,t)] 



mo)) 



where we have used the usual time-ordering operator 

t H. 



the quadratic Hamiltonian accounting for interactions 
effects. 



III. OPTICAL PROPAGATION OF MATTER 
WAVES: THE ABCD THEOREM. 



The strategy exposed in this paper consists in pick- 
ing up, among the possible polynomials P, the element 
which minimizes an appropriate distance measure to the 
mean-field potential. In geometric terms, this polynomial 
appears as the projection of the mean-field potential onto 
the vector space spanned by second-order polynomials in 
position and momentum. This potential will be referred 
to as the "lensing potential" , denomination which will be 
justified in Sec. HVl We define a distance analogous to 
the error function used in (3fT |. which involves the poly- 
nomial P and the quantum state | </>(£)) resulting from the 



This section gives a remainder on a general result 
- called the ABCD theorem - concerning the propa- 
gation of matter waves in a time-dependent quadratic 
potential, which is the atomic counterpart of the ray 
matrix formalism frequently used in optics [2(| . It shows 
that the evolution of a Gaussian wave-function under 
an Hamiltonian quadratic in position and momentum 
is similar to the propagation of a Gaussian mode of 
the electric field in a linear optical system. A detailed 
description of this theoretical result of atom optics is 
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given in the references [25|, l35| . 

One considers a time-dependent quadratic Hamilto- 
nian such as: 



H a + Pi(r,p,t) 



Pf3(t)f> 
2m 



ipa(t)f - ~rS(t)p 



-ry(t)r - mg(t) ■ f + f (t) • p + h(t) 



(6) 



a(t), (3(t), -f(t) and 5(t) are 3x3 matrices [ii] ; f(t) and 
g(i) are three-dimensional vectors; h(t) is a scalar, and~ 
stands for the transposition. Here we use this Hamilto- 
nian to approximate the nonlinear Hamiltonian ([1]) . The 
Hamiltonian ^ is indeed appropriate to describe several 
physical effects [H, [37j ■ 



ABCD propagation of a Gaussian 
wave- function. 



The propagation of a Gaussian wave-packet in such 
an Hamiltonian can be described simply as follows. Let 
<p(r, t) be an atomic wave packet initially given by: 



(r, to 



VldctXo 



■e 2h 



(r-r c0 )Y X a (r— r o0 )+-f PcO'(r— r c o) 



(7) 



The 3x3 complex matrices Xq, Yq represent the ini- 
tial width of the wave packet in position and momen- 
tum respectively: X = iD(Ax(t ), Ay(t ), Az(t )) and 
Y = D(Ap x (t Q ),Apy(t Q ),Ap z (to)), with D standing for 
a diagonal matrix. The vectors r c o, p c o give the initial 
average position and momentum. The ABCD theorem 
for matter waves states that, at any time t > to, the 
wave-packet (j>(r, t) satisfies: 

pjrS{t, t , r e0 ,p c0 ) , 
<h( t t ) = . p f(r-r rf )y f X t - 1 (r-r ct ) + lp cl .(r-r cl ) 

S(t, to, r c o, Pco) is the classical action evaluated between 
t and to of a point-like particle which motion follows the 
classical Hamiltonian H (r, p, t) and with respective ini- 
tial position and momentum r c o, p c o- The width matrices 
in position X t and momentum Y t , and the average po- 
sition and momentum r ct , p ct at time t are determined 
through the same 6x6 "ABCD" matrix: 



r c t 

x, 

Yi 



A{t,to) B{t,t ) 

C{t,to) D(t,t ) 

A(t,t Q ) B(t,t ) 

C(t,t ) D(t,to) 



-PcO 

X 
Y Q 



( S(t,to) 



The ABCD matrix -noted compactly M(t. tn)- and the 
vectors ^, <j> can be expressed formally as [371 ]: 



M(t,t ) 



<j>(t,t ) 



T 



exp 



, f ,(a(t') f3(t>) 
1 lit') S(t') 



d£ M(t,t') 



f(f) 
g(i') 



(8) 
(9) 



Although the former expressions seem rather involved, in 
all cases of practical interest, the ABCD^cj) parameters 
can be determined analytically or at least by efficient 
numerical methods. 



B. Interpretation of the ABCD propagation and 
aberrationless approximation. 

The phase-space propagation provides a relevant in- 
sight in the transformation operated by the ABCD ma- 
trix. Consider the Wigner distribution of a single-particle 
density operator evolving under the Hamiltonian (|6j). 
The Wigner distribution at time t is related to the dis- 
tribution at time £q by the following map: 



W(r,p,t) = W( D(v-C)--B(p- 

\ 771 



nuj>) , 



[ii] S(t) = —a(t) to ensure the Hamiltonian hcrmiticity 



-mC(r - + ^(P - m <i>) > to 



where the matrices A, B, C, D and vectors £, <j> are again 
evaluated at the couple of instants (t, to). The action of 
the evolution operator onto the Wigner distribution is 
thus amenable to a time-dependent linear map. The fact 
that ABCD matrices are symplectic [25| implies that 
this map is unitary: such evolution preserves the global 
phase-space volume, and the quality factor of an atomic 
beam in the sense of (38j . 

In photon as in atom optics, the aberrationless ap- 
proximation consists in assuming that the Gaussian func- 
tion ([7]) is a self-similar solution of the propagation equa- 
tion in spite of the non-linearity, the evolution of which 
is given by the ABCD propagation. The propagation 
is thus described through a map which preserves the 
phase-space density. This is an approximation, since 
for atomic or light beams evolving in nonlinear media, 
the phase-space density indeed changes during the prop- 
agation. Nonetheless, the aberration- less approximation 
is reasonable for sufficiently diluted clouds, subject to a 
weak mean-field interaction term, for which an initially 
Gaussian wave-function will not couple significantly to 
higher-order modes. Furthermore, this approximation in 
atom optics is entirely analogous to the aberration-free 
treatment realized in non- linear optics, the predictions 
of which concerning the width evolution of a light beam 
have been verified experimentally [391 ] . One can thus ex- 
pect that the aberrationless approximation will consti- 
tute a good description of the propagation in atom optics 
as well. Indeed, the validity of the aberrationless approx- 
imation will be confirmed in Sec. IV Di on the example of 
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a gravitational atomic resonator: its predictions on the 
sample size evolution are in good agreement with those 
of a paraxial treatment of the wave-function propagation 
which does not assume the preservation of a Gaussian 
shape. 



IV. ABCD MATRIX OF A FREE-FALLING 
INTERACTING ATOMIC CLOUD. 

Let us apply the method discussed above to describe 
the propagation of a free-falling Gaussian atomic wave- 
packet. In the aberrationless approximation, such a 
wave-packet is simply determined by the parameters 
ABCD^cf) and by the phase associated with the action. 
In view of the resonator stability analysis, we will focus 
on the computation of the ABCD matrix in presence 
of the mean-field potential. We consider only the 
leading-order nonlinear corrections, associated with the 
first-order lensing polynomial P, (r, p, t). 

This section begins with the determination of this po- 
tential defined by Eq.©. A formal expression of the 
atom-optical ABCD matrix, taking into account this 
lensing potential, is obtained. An infinitesimal expan- 
sion of this expression shows that the mean-field inter- 
actions effectively play the role of a divergent lens: the 
atom-optical ABCD matrix of the free-falling cloud evo- 
lution is similar to the optical ABCD matrix associated 
with the propagation of a light ray through a series of in- 
finitesimal divergent lenses. In our case, the propagation 
axis is the time, and the infinitesimal lenses correspond 
to the action of the mean-field potential in infinitesimal 
time slices. 



A. Determination of the lensing potential. 



The square of the free-evolving wave- function thus reads: 



l0 (O) (r,t)| 2 



-3/2 



W xt WytW zt 



(*-*ctr (v-vctr (z-zctr 
— 

_ g xt yt zt 



We use this expression to determine the first-order lens- 
ing polynomial P/^(r,p,i). Since this operator acts on 
Gaussian wave-functions, differentiation is equivalent to 
the multiplication by a position coordinate, so the action 
of the momentum operator is indeed equivalent to that 
of the position operator up to a multiplicative constant. 
One can thus, without any loss of generality, search for 
a lensing polynomial P, (r, t) involving only the posi- 
tion operator. With this choice, the error function ([3|) 
minimized by the polynomial P becomes simply: 

P(P(t),|^ 0) (t)}) = J d s r|0 (o) (r,t)| 2 (P(r,t)~ gi \t m ( r ,t)\- 

Expanding the polynomial P/ 1 '' (r, t) around the central 
position r ct , a parity argument shows that the linear 
terms vanish: 

P ; (1) (r,i) = 9I [c (t) - c x (t)(x - x ct ) 2 

- c y(t)(y ~ Vet) 2 ~ c z (t)(z - z ct ) 2 ] 

By definition of the lensing polynomial, the error function 
(fl~2l must be stationary with respect to each coefficient 
Cx,y,z,o(t) ! which leads to: 



co(t) 
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with V(t) = {2i:f/ 2 w xt w v tw zt (12) 



Only the quadratic term intervene in the ABCD matrix: 
the coefficient co(t) merely adds a global additional phase 
to the wave-function, which does not change the subse- 
quent stability analysis. 



We assume that the condensate, evolving in the Hamil- 
tonian {1} , is initially at rest and described by a Gaussian 
wave- function: 



(x,y,z,t ) 



r-3/4 



(10) 



y/W xQ WyQW z0 

It is easy to show that, when one considers the 
interaction-free evolution, the widths are given at time 
t > t by: 



n 2 



Wit 



-(*-*>)* 



(11) 



'70 



for i — x,y, z. This result can be easily re- 
trieved by considering the initial width matrices Xq — 
iD(w x0 ,w y0 ,w z0 ) and Y = -^T)(l/w x0 ,l/w y0 ,l/w z0 ) 
for the wave- function, and applying the free ABCD ma- 
trix H: 



A(t,t ) B(t,t ) 
C(t,t ) D(t,to) 



1 t-t 
1 



B. Formal expression of the effective ABCD 
matrix. 

We can readily express the ABCD matrix associated 
with the evolution under H^(t). Writing this Hamilto- 
nian in the form of Eq. ([S]), and using the formal expres- 
sion ([S]) of the ABCD matrix as a time-ordered series, 
one obtains: 



MW(t,t ,X ) 



T 



exp 



dt' 



a(f) (3(f) 
7 (t') 8(f) 



(13) 

In contrast to the usual linear ABCD matrices, this ma- 
trix now depends on the input vector through the initial 
position width matrix Xq [iii] . A brief inspection of Eq. 



[iii] The Hamiltonian H^(t) and the lensing polynomial P,^(r,i) 
depend of course also on Xo, but we do not mention this depen- 
dence explicitly to alleviate the notations. 
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© and of the Hamiltonian i?W (t) 
H (1) (t) = ^- + mgz + gi[c {t) 



c y (t)(y - yet? 



c x {t)(x - 
c z (t){z - z ct ) 2 ] 



Xct 



(14) 

shows that the matrices in the exponential read a(t) = 
S(t) = 0, (3(t) = 1 and 7 (i) = ^D^^^^.c,^)). 
Using Eq. (TH?)) . one readily obtains the elements of the 
quadratic matrix 7: 

1 



7«(*) = 



5j 



m wf t (w xt w vt w zt ) 



(15) 



for i = x,y, z, with the widths w x , y ,zt given by Eq. (fTTj) . 
A significant simplification arises because j(t) is diago- 
nal: one needs only to compute the exponential of three 
2x2 matrices associated with the orthogonal directions 
O x ,O y ,O z . The ABCD matrix is simply the tensor prod- 
uct of those: 



MM(t,to,X )=®i= w , y , z T 



exp 



dt' 



la 





7«(*) 



(16) 



C. Propagation in a a series of infinitesimal lenses. 



An infinitesimal expansion of (| 18[) shows that the evo- 
lution between t and t + dt is described by the ABCD 
matrix: 

1 dt 



M^(t + dt,t,X ) 



(17) 



j(t)dt 1 

It can be rewritten as a product of two ABCD matrices: 



M (1 \t + dt,t,X) 



1 dt 
1 



1 







j(t)dt 1 



(18) 



If these were 2x2 matrices, in the optical formalism, the 
first matrix would be associated with the propagation of 
a ray on the length dt and the second matrix, of the form 



1 

-dt/f 1 



(19) 



would model a lens of infinitesimal curvature dt/f. One 
can thus consider, by analogy, that this second 6x6 ma- 
trix realizes an atom-optical lens which curvature is the 
infinitesimal 3x3 matrix T)("f xx (t), Jyy(t), J zz {t)) dt. Be- 
sides, one can exploit the fact that it is a tensor product: 
if one considers each direction O x ,O y , O z separately, the 
propagation amounts - as in optics - to a product of 2 x 2 
matrices, which makes the analogy with a lens even more 
transparent. The resulting 6x6 ABCD matrix is simply 
given by the tensor product of those. Transverse de- 
grees of freedom are, nonetheless, coupled to each other 
through the lensing potential. It is worth noticing that 
the focal lengths f x , f y , f z have here the dimension of a 
time, and are negative if one considers repulsive interac- 
tions: the quadratic potential Pf\r, p, t) acts as a series 
of diverging lenses associated with each infinitesimal time 
slice. 



D. Expression of the nonlinear ABCD matrix with 
the Magnus Expansion. 

Because of the time-dependence of the Hamiltonian 
H^'(t), the time-ordered exponential in (fig)) cannot, in 
general, be expressed analytically. Fortunately, a useful 
expression is provided by the Magnus expansion [40| : 



M^(t,t ,X ) = 



^i—x : y,z 



exp 



la 



, dh / dt 2 [N i (t 1 ),N i (t 2 )} 
*Jo Jo 



with Ni(t) 



1 

m(t) 



(20) 



This 



where ^>i— Xt y,z denotes again a tensor product, 
expansion has the advantage to preserve the unitarity 
of the evolution operator: at any order, the operator 
obtained by truncating the series in the exponential is 
unitary. The Magnus expansion can be considered as 
the continuous generalization of the Baker-Hausdorff 
formula [4l| giving the exponential of a sum of two op- 
erators A and B as a function of a series of commutators 
along exp(A + B) = exp AexpBexp([A, B}/2).... The 
Magnus expansion has been successfully applied to solve 
various physical problems, among which differential 
equations in classical and quantum mechanics [42l |. 
spectral line broadening [43| . nuclear magnetic reso- 
nance [441, multiple photon absorption |45ll and strong 
field effects in saturation spectroscopy |46j |. 

The first-order term £li(t,to) in the argument of the 
exponential can be expressed as 



ni(t,to) = jT dt'N(t')= ( 





( 7 )r 



(21) 



with the duration r = t— to and the average quadratic di- 
agonal matrix (7)^ = 1/r /. dtju(t). Exact expressions 
for {j)u are given in Eq. (|C1[) of Appendix [Cl for a cylin- 
drical condensate. Without this symmetry, the matrix 
elements (7)^ cannot be evaluated analytically to our 
knowledge, but are nonetheless accessible with efficient 
numerical methods [iv] . The first-order ABCD matrix 

M^fafaXo) = reads [v] : 



M 



(i) 



cosh(( 7 ) 1 / 2 T) (7)~ 1/2 sinh((7) 1 / 2 r) 
(7) 1/2 sinh(( 7 ) 1 /2 r ) cosh(( 7 ) 1 /2 r ) 



(22) 



[iv] In the short expansion limit considered later where |t — £0 1 <S 
mwf(to)/h, the average quantities {yu) can be approximated by 
the instantaneous value of the quadratic coefficient 7^ at the 
center of the considered time interval. 

[v] For repulsive interactions, all the eigenvalues of the matrix 7 
are positive, and by convention its square root has also positive 
eigenvalues. 
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As expected, this main contribution of the Magnus 
expansion is independent of the ordering of the suc- 
cessive infinitesimal lenses, and can be interpreted as 
the ABCD matrix of a thick lens with finite curvature. 
This expression is similar to the paraxial ABCD matrix 
obtained in [l8| to describe the interactions between an 
atom laser and a condensate of known wave-function. 

In the following developments, we use mainly this first- 
order contribution to the Magnus expansion. In order to 
justify this approximation, we have performed a second- 
order computation of the ABCD matrix M^- 1 ' (t, to, Xq) 
in Appendix iBl This second-order correction is weighted 
by the small parameter e = (r/r c ) 4 , depending on the 
ratio of the duration r = t — to to a time-scale t c , which 
reads for a spherical cloud of radius Wq: 



T < T C 



/ Wo \ x / 6 mwq 



\47ra/ 



(23) 



One checks that the first-order expansion is valid for 
an arbitrary long time (r c — > oo) as interaction effects 
vanish (a — > 0). Considering a sample of initial radius 
wq = 10 /im, and using the s-wave scattering length 
a ~ 5.7nm of the 87 Rb [Hj], one obtains r c = 0.31 s. 
The convergence of the Magnus series is indeed guaran- 
teed when the following inequality is satisfied [4l| : 



N„ 



6£ \\N(t')\\ < ln(2), 



(24) 



to 



and our second-order computation gives an additional 
heuristic indication of convergence for a flight duration 

t < T c . 



V. STABILITY ANALYSIS OF A 
MATTER- WAVE RESONATOR 

In this Section, we apply the method of the ABCD 
matrix to discuss the propagation of an atomic sample 
with mean-field repulsive interactions in a matter-wave 
resonator [471 ]. The considered resonator involves a se- 
ries of focusing atomic mirrors. In this system, there is a 
competition between the transverse sample confinement 
provided by the mirrors and the expansion induced by 
the repulsive interactions, which determines the maxi- 
mum size of the sample during its propagation. In order 
to keep the sample within the resonator, its transverse 
size must stay smaller than the diameter of the laser 
beams realizing the atomic mirrors. If this criterium 
is met during the successive bounces, the resonator is 
considered as stable. The ABCD matrix method devel- 
oped previously, giving an easy derivation of the sample 
width evolution, is well-suited to discuss this issue. One 
assumes an initial Gaussian profile for the sample wave- 
function. The atomic wave propagation in-between the 
mirrors is treated in the aberrationless approximation, 
and described by the nonlinear ABCD matrix (|2"2"|) ac- 
counting for self-interaction effects. The evolution of the 



sample width obtained with this method is compared to 
the behavior expected from a non-perturbative paraxial 
approach. 



A. Resonator description 

The considered matter- wave resonator is based on the 
levitation of a free-falling two-level atomic sample by pe- 
riodic vertical Raman light pulses. This proposal is de- 
scribed in detail in the reference , but we remind here 
its main features for the sake of clarity. In the absence of 
light field, the atomic sample propagates in the Hamil- 
tonian (TT]). We consider an elementary sequence which 
consists in a pair of two successive short vertical Raman 
7r pulses [H. Each pulse is performed by two counter- 
propagating laser beams of respective frequencies uj up , 
Udown and wave- vectors k up = kz, kdown = —kz equal 
in norm to a very good approximation and of opposite 
orientation. The first Raman pulse propagates upward 
with an effective vertical wave-vector k e ,i = 2kz and 
corresponds to laser frequencies uj up — u>2, ujdown = <^i', 
the second one propagates downward with an effective 
vertical wave-vector k e> 2 = —2kz and with the laser fre- 
quencies LO U p — UJ3 and uidown = w 4- The frequencies 
^1.2,3.4 are adjusted so that both Raman pulses have the 
same effective frequency u> e = \oj up — LOdown\, satisfying 
the resonance condition (47j u> e = UJ2 — u>i = 0J4 — 0J3 = 
LUb a — 2hk 2 /(mh). The intermediate level involved dur- 
ing the Raman pulses (of energy E — h(uj a + u> e )) is 
taken sufficiently far-detuned from the other atomic en- 
ergy levels to make spontaneous emission negligible [vi] 
. After adiabatic elimination of the intermediate level, 
the action of the Raman pulses can be modelled by the 
effective dipolar Hamiltonian: 



H d ip(t) = -Klba(T,t) COs(u e t - k e l,2 ' ? ) (|6)(o| + 



a)(b\) 
. (25) 

Each pair of pulses acts as an atomic mirror, bringing 
back the atoms in their initial internal state a, and pro- 
viding them with a net momentum transfer of Ap ~ 4/ik. 
The atomic motion is sketched on Fig. [T] in the energy- 
momentum picture. 

This sequence can be repeated many times. If the period 
T in-between two successive atomic mirrors is set to 



T := Tn = 



Ahk 
mg 



(26) 



the acceleration provided by the Raman pulses compen- 
sates on average that of gravity: the cloud levitates 
and evolves inside a matter- wave resonator [47[. An 
analogous system has been realized experimentally re- 
cently m. 



[vi] In practice, a detuning on the order of the GHz - experimentally 
compatible with 7r pulse of duration shorter than the ms [48j - is 
sufficient to discard spontaneous emission. 
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|z=0,b> 



Total Energy 
E = f(p) 




|z=z n ,a> YPi 



|z=0,a > 



Momentum p 



-2hk 2nk 



FIG. 1: (Color online) Evolution of the atomic sample in the 
energy-momentum picture. The total energy includes the ki- 
netic, gravitational and internal energy. The atoms are ini- 
tially at rest (p = 0), at the altitude zo, and in the lower 
state a. The starting point is thus at the intersection of the 
paraboloid (a,zo) and of the vertical axis (p=0). In between 
the pulses, the motion of the atomic sample in the gravity 
field is conservative: it corresponds to a leftward horizontal 
trajectory of the representative point. 



B. Focusing with atomic mirrors. 

Matter-wave focusing can be obtained, in principle, 
with laser waves of quadratic intensity profile [HoL l5l| 
or alternatively of spherical wave- front [47| . We concen- 
trate on the focusing obtained with an electric field of 
quadratic intensity profile [HTj] , the discussion of which 
is less technical. The Rabi frequency considered for the 
Raman pulses of the resonator depends quadratically on 
the distance to the propagation axis O z [vii] : 



1 



V 



2wi 



n (t) (27) 



These Raman pulses generate a quadratic position- 
dependent light-shift proportional to the field intensity 
and thus to the square of the Rabi frequency (f2"T)) . Af- 
ter the pulse, the atomic wave-function initially in of the 
form of Eq. ([7]) is thus multiplied by a factor yielding the 
input-output relation: 



ipout(r,t) 



— g«2fc(2 



-Zo) e -i(x +y )/w u 



»e^Vin(r,t) (28) 



with 4>' Q a constant phase added at the condensate cen- 
ter r = (0, 0, Zo) during the pulse. The outgoing wave- 
function can thus be put again in the form of Eq. (|7|) if 
one replaces po by p\ = po + 2hk, and Xq, Yq with 



Xi 



D(- 



h 

-!//,-!//, 0) 



A 
Y 



(29) 



[vii] Close to the propagation axis, this quadratic profile can be repro- 
duced to a good approximation with Raman pulses of Gaussian 
intensity profile. 



/3,03 are the 3x3 identity matrix and null matrix, 
D(—l/f, —1/f, 0) is as previously a 3x3 diagonal matrix. 
The focal time is: 



/ = 



u las 



2h 



(30) 



Eq. shows that the pulse acts as a lens in the 



transverse directions O x , O y [viii] 



The strength of the focusing which can be achieved 
with such atomic mirrors [ix] is indeed limited by the 
quasi-uniformity required for the Rabi frequency on the 
condensate surface, in order to perform an efficient pop- 
ulation transfer with the Raman 7r-pulse. Considering a 
cigar-shaped cloud of small width w± along the O x , O y 
axis, one may require that the Rabi frequency difference 
between the border and the center of the cloud satisfies: 
\fl(w±,0,z,i) - n (t)\/\n (t)\ < e. This yields readily a 
lower bound on the focal time /: 



/ > 



2he 



(31) 



With a reasonable bound of e = 10 -2 , a cylindrical cloud 
of 87 Rb atoms of transverse size w± ~ 10/im, one obtains 
a minimum focusing time: / > 6.7s. A back-on-the- 
envelope computation of the reflection coefficient shows 
that the losses resulting from such an inhomogeneity of 
the Rabi frequency are on the order of 10~ 3 . 



C. Resonator stability analysis. 

We now investigate the non-linear ABCD propagation 
of a cigar-shaped sample in the resonator. As a specific 
example, we consider a cloud of 87 Rb atoms taken 
in the two internal levels \a) = \5Si/2,F — 1) and 
\b) = \5Si/2,F = 2). In-between the Raman mirrors, 
the whole sample is expected to propagate in the 
ground state \a). We consider a sample of N = 10 5 
atoms, of initial dimensions w x = w y = w r = 10 \xm 
and w z — 100 /im, and we use the s-wave scattering 
length a ~ 5.7 nm of the Rubidium. We investigate 
the evolution of this sample during a thousand bounces 
and for various mirror focal times. Keeping a significant 
atomic population inside a matter-wave resonator during 
such a big number of reflections is challenging, but not 
impossible in principle given the high population transfer 
which has been achieved experimentally with Raman 



[viii] The absence of focusing in the direction of laser beam propaga- 
tion O z is not critical since it does not drive the cloud out of the 
beam. 

[ix] The considered atomic mirrors consist indeed not in a single, but 
in a double Raman pulse. This does not change the qualitative 
discussion of this paragraph. 
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pulses [52j [x] One obtains the value To ~ 1.5 ms 
for the period between the Raman mirrors. This time 
scale is much shorter than the duration r c ~ 0.3 s 
found for the validity of the first-order Magnus ex- 
pansion associated with a spherical cloud of radius 
wo = 10 /im. This shows that the ABCD matrix of 
the cigar-shaped condensate is well-approximated by 
the leading order [Eq. ([22)1 ] of the Magnus expan- 
sion [xi] . Furthermore, the free-propagation time To 
is also much shorter than the time-scale r r = mw^/h 
associated with the free expansion of the transverse 
width, so that one can safely approximate the average 
quadratic coefficient (7) with the instantaneous value 

(7) ~ l(w xTo /2,W y T /2,W zTo/2 ). 

To compute the evolution of the transverse and lon- 
gitudinal sample width, one proceeds as follows. As in 
Section [IV A[ one starts with initial width matrices Xq = 
iD(w x0 ,w y o,w z o) and Y Q = ^B(l/w x0 ,l/w y Q,l/w z0 ) 
and computes the interacting ABCD matrix (|22|) as a 
function of these initial widths. During the first cycle, 
one multiplies the corresponding vector (Xo,Yq) succes- 
sively with nonlinear ABCD matrix (|22|) and with the 
mirror ABCD matrix (|29|l . The new width matrices 
(Xi, Yi) are obtained, from which one can infer the non- 
linear ABCD matrix for the next propagation stage. The 
iteration of these algebraic operations is a straightfor- 
ward numerical task. The results, depicted on Fig. [21 
show that the transverse width oscillates with an ampli- 
tude and a period which both increase with the mirror 
focal time. The maximum sample size is w r = 25 /im 
and w r — 60 /im for the respective focal times / = 20 s 
and / = 100 s. Considering for instance a laser beam of 
waist w = 100 fim in the experiment, one sees that with 
those focal times the atomic cloud remains within the 
light beam and is thus efficiently confined transversally 
in the resonator. As expected, the use of Raman mirrors 
with a stronger curvature allows one to shrink the trans- 
verse size of the stabilized cloud. Fig. [3] shows the evolu- 
tion of the maximum sample transverse size as a function 
of the mirror focal time. The extended ABCD matrix 



analysis presented in this paper allows thus to determine 
efficiently the minimum amount of focusing required to 
keep the sample within the diameter of the considered 
Raman lasers. In that respect it can be used to optimize 
the trade-off, exposed in the previous paragraph, between 
strongly focusing or highly reflecting atomic mirrors. 

Width (pm) 

200l 1 1 1 1 1 1 1 1 — ^ 1 




°0 100 200 300 400 500 600 700 800 900 1000 
Number of bounces 



FIG. 2: (Color online) Evolution of the transverse and longitu- 
dinal width of the sample (/im) during the successive bounces 
in the cavity (numbered from 1 to 1000), for the mirror focal 
times / = 20 s (blue), / = 50 s (green) and / = 100 s (red). 
The dashed line represents the evolution of the transverse 
width in the absence of focusing with the Raman mirrors. 
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FIG. 3: (Color online) Maximum sample transverse width 
(pm) during the evolution in the resonator as a function of 
the Raman mirror focal time (s). We have considered the first 
1000 bounces to determine this maximum. 



[x] We treat the wave-propagation in the resonator as if the atomic 
cloud was entirely reflected on the successive atomic mirrors. In- 
deed, even if resonant Raman pulses can perform a population 
transfer with an efficiency close to 99% lo2H . the residual losses 
become significant after a big number of bounces in a real ex- 
periment. This results in a gradual decrease of the mean-field 
interactions, which could be accounted for in a more sophisti- 
cated model. Our point here is simply to illustrate the nonlinear 
ABCD method on a thought experiment, and we thus adopted 
a simplified approach with perfect atomic mirrors. 

[xi] We have computed the time-scale t c determining the validity 
of the first-order Magnus term for spherical wave-packets only. 
Nonetheless, a basic dimensional analysis shows that for a cigar- 
shaped cloud, the time-scale determining the validity of the first- 
order Magnus term is bounded below by the time T c given by 
Eq. (1236 and computed by setting wq equal to the smallest cigar 
dimension. 



D. Comparison with the predictions of the 
nonlinear paraxial equation. 

As exposed in Appendix [Al the propagation of an 
atomic beam with a longitudinal momentum much 
greater than the transverse momenta can be alter- 
natively described by a paraxial wave equation of 
the form (|A2[) . Furthermore, if the linear density of 
the atomic beam is uniform, the nonlinear coefficient 
intervening in this paraxial equation is a constant. As 
in nonlinear optics [53| and in 2D condensates [54| . 
this equation induces a universal behavior in paraxial 
atomic beams [Hj]: the transverse width oscillates 
with a frequency independent from the strength of the 
interaction. The width oscillations, depicted on Fig. [21 
indeed allow one to confront the results of our method, 



9 



which uses a non-paraxial wave equation treated in the 
aberrationless approximation, to the predictions of the 
full nonlinear paraxial equation with a uniform nonlinear 
coefficient. We stress that this second approach leaves 
the nonlinear term as such and does not assume that 
the Gaussian shape of the atomic beam is preserved. In 
this sense it is more exact than the radius of curvature 
method used in Appendix [XJ It is also approximate, 
since the atomic beam is neither paraxial nor of uniform 
linear density. Nevertheless, it is remarkable that both 
treatments agree very well on the oscillation period of 
the width. 

To apply the paraxial description, one models the ac- 
tion of the successive mirrors on the transverse wave- 
function with an average potential. The lens operated 
by each Raman mirror, of focal time /, imprints a phase 
factor of e l W r [see Eq. ([28]) and Eq. ([30])]. The series 
of lenses, separated by the duration To, thus mimics the 
following effective quadratic potential: 



V 



-L,lens 



2h 2 T f 



(32) 



Let us consider the nonlinear contribution, given 
by a contact term of the form V±.i n t(r) = 
g I \ip//(z)\ 2 \ip ± (r)\ 2 ip ± (r), with ij)//(z) the longitudi- 
nal wave-function [Eq. (|A1[) of Appendix [A] . The term 
gi\ip//(z)\ 2 appears as an effective nonlinear coupling 
coefficient for the transverse wave-function depending 
on the altitude z. Adding this nonlinear contribution 
to Eq. (|Aip . one obtains a 2D nonlinear Schrodinger 
equation (NLSE): 



ihd ( i)j_(x,y,C) = 



fr 2 

-f-(flg 



- d 2 ) + gA^OV^i 



2h 2 T Q f 



ip±(x,y,C), (33) 



£ is a parameter defined in Eq. (|Al[) equivalent to the 
propagation time, a the scattering length and r 2 = x 2 + 
y 2 . Setting K = one can recast this equation in the 
same form as in [53] where the propagation of a light wave 
in a quadratic graded index medium was considered: 



2iKd c i/j A 



-d 2 



(87ra|%| 2 ) |V> 



K 2 



1 

Wo 



We now make the assumption that the variations of the 
non-linear coefficient 8ira\ip//(£)\ 2 with £ are sufficiently 
smooth to have a negligible impact on the period of the 
sample width oscillations. This assumption seems rea- 
sonable for the considered cigar-shaped cloud, which has 
a slow longitudinal expansion in comparison with the os- 
cillation period [see Fig. [2]. This hypothesis is indeed 
validated a posteriori, since it leads to predictions in ex- 
cellent agreement with the results of the ABCD method 
discussed above. Once the nonlinear coefficient is approx- 
imated with a constant, one can readily apply the results 



derived in [53j, |55j , which show that Eq. ([34)1 yields trans- 
verse oscillations of universal frequency: 



(35) 



The results obtained from the perturbative ABCD ap- 
proach are confronted with this prediction on Fig. 3] The 
agreement improves as the mirror focal time increases, 
and it is in fact already good (4%) for a focal time of 
/ = 3 s and attains 0.7% for a focal time of / = 50 s. 
As discussed above, focal times shorter than / = 20 s 
seem incompatible with the reflection coefficient desired 
for the atomic mirrors. The disagreement observed be- 
low / < 3 s may be attributed to a failure of the paraxial 
approximation to describe the propagation of the sample 
in our system. 




10 15 20 25 30 35 40 45 50 
Mirror focal time (s) 



FIG. 4: (Color online) Period of the transverse width oscil- 
lations (s) in the matter-wave resonator as a function of the 
Raman mirror focal time (s). The full and the dashed line give 
the oscillation periods obtained respectively through the per- 
turbative ABCD approach and through the nonlinear paraxial 
wave equation. 



VI. CONCLUSION 

This paper exposed a treatment of the non-linear 
Schrodinger equation involving theoretical tools from op- 
tics and atom-optics. The ABCD propagation method 
for matter waves has been extended beyond the linear 
regime thanks to a perturbative analysis relying on an 
atom-optical aberrationless approximation. We have de- 
rived approximate analytical expressions for the ABCD 
matrix of an interacting atomic cloud thanks to a Mag- 
nus expansion. This matrix analysis has been applied to 
discuss the propagation of an atomic sample in a perfect 
matter- wave resonator. We have shown that such sam- 
ple can be efficiently stabilized thanks to focusing atomic 
mirrors. We have found that the nonlinear ABCD prop- 
agation reproduces to a good level of accuracy the uni- 
versal oscillations expected from the nonlinear paraxial 
equation for matter waves (55| , which makes it a promis- 
ing tool to model future nonlinear atom optics exper- 
iments and a seducing alternative to previous numeri- 
cal methods applied to matter- wave resonators [50(. We 



10 



have also highlighted an other optical method, involv- 
ing more stringent assumptions - paraxial propagation, 
cylindrical symmetry and constant longitudinal velocity 
- and also relying on the aberrationless approximation. 
This last method enables one to address self-interaction 
effects in the free propagation through a complex param- 
eter [defined in Eq. (|A13[) ]. which is analogous to a radius 
of curvature, and the evolution of which is very simple 
[Eq. (|A14[) ]. As far as the beam width is concerned, the 
effect of self-interactions can be interpreted as a scal- 
ing transformation of the free propagation by a factor 
depending on the matter- wave flux T [See Eq. (|A15|) ]. 
Both approaches are relevant to study interaction effects 
on the stability of atomic sensors resting on Bloch oscil- 
lations [56j ] , on the sample propagation in coherent inter- 
ferometers [57j . An interesting continuation of this work 
would be to develop a nonlinear ABCD matrix analysis 
beyond the aberrationless approximation. 
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APPENDIX A: THE METHOD OF THE 
NON-LINEAR RADIUS OF CURVATURE. 

This method adresses the paraxial propagation of a 
monochromatic and cylindrical matter-wave beam. It 
relies on the introduction of an effective complex radius 
of curvature [23, 26], which evolution is especially simple, 
even for a self-interacting beam. It has been applied 
successfully by Belanger and Pare [24[ to describe 
self focusing phenomena of cylindrical optical beams 
propagating in the paraxial approximation, and it works 
equally well for matter waves propagating in the same 
regime. This is typically the case for an atom laser beam 
falling into the gravity field, for which the transverse 
momentum components become negligible compared to 
the vertical momentum after sufficient time 181. 



We consider a mono-energetic wave-packet propagat- 
ing in the paraxial regime, and evolving in the sum 
of a longitudinal potential V//(z) and a transverse one 
V±(x, y, z), which may also vary slowly with the longitu- 
dinal coordinate z. This section begins with a brief re- 
mainder on the paraxial equation for matter waves 33 1 , 
and on its spherical- wave solutions in the linear case 23 1 . 



the nonlinear propagation [23j | , at the cost of certain ap- 
proximations, and thanks to the introduction of a gen- 
eralized radius of curvature depending on the coupling 
strength. Our treatement of the nonlinear matter wave 
propagation follows step by step the approach of Belanger 
and Pare for optical waves [24| . 



1. The Paraxial Equation for Matter Waves. 

Our derivation of the nonlinear paraxial wave-equation 
follows the treatment done in [33| . The wave- function is 
factorized into a transverse and longitudinal component: 

2p(x,y,z) = ip ± (x,y,z)ip // (z) 

The longitudinal component obeys a ID time- 
independent Schrodinger equation, 



h 2 



2m dz 2 

which can be solved with the WKB method: 



%j}f/{z) 



p(z) 



exp 



dup(u) 



T = J d 2 rj_£^|i/)(r_L, z)\ 2 is the atomic flux evalu- 
ated through any infinite transverse plane, the trans- 
verse wave-function ip±_ being normalized to unity 

J d 2 r±\i>±{r ± , z)\ 2 = 1. p(z) = ^2m(E-V/,{z)) is the 
classical momentum along z, and z is the associated clas- 
sical turning point verifying p(z ) = 0. The transverse 
wave-function i/jx, assumed to depend slowly enough on 
the coordinate z to make its second derivative negligible, 
verifies the equation: 



in a. 



2m 



9y) 



Vj_(x,y,z) 



ipj_(x,y,z) = 0. 



This equation can be simplified with a variable change 
in which the longitudinal coordinate z is replaced by the 
parameter £: 



COO 



(Al) 



which corresponds to the time needed classically to prop- 
agate from the turning point zq to the coordinate z [xii] 
. The wave equation becomes 



ififif + f- (d 2 x + d 2 y )-V ± (x,y,0 
2m 



^±{x,y,0 = 0, 



(A2) 



It is remarkable that such solutions can be extended to 



[xii] Indeed, this parameter appear as proportional to the proper time 
experienced by the atom on the classical trajectory determined 

bypO)[H. 
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We assume from now on that the transverse potential 
V±(x,y, z) has a cylindrical symmetry. If one sets K — 
m/h and V±(x,y,Q = ^K 2 {Qr 2 with r 2 = x 2 + y 2 , 
Eq. (|A2[) has the same form as the paraxial equation for 
the electric field used in [241 ]: 

[3 2 T + 2iKd ( - KK 2 {Qr 2 } ^(x, y, C) = , (A3) 

It is worth noticing that, as a consequence of our variable 
change, the derivative with respect to the longitudinal 
coordinate z has been replaced by a time derivative with 
respect to C- 



2. Spherical Wave solutions to the linear equation. 



One looks for solutions of Eq. (|A3j) of the kind: 

K 



i>±{x,y,C) = -4(0 exp 



25(0 



(A4) 



with again K = m/h. Such function is a solution if 
and only if the parameter g(£) - called complex radius of 
curvature, and homogenous to a time for matter waves - 
satisfies the following equation: 



q 2 K 
and if the amplitude A(() verifies: 
A' 1 











(A5) 



(A6) 



The prime stands for the derivative with respect to 
(,. In the absence of the transverse potential, i.e. 
V±(x,y, C) = 0, an obvious evolution is obtained with 

?(0 = C 

These equations imply a relation between the ampli- 
tude and width of the wave- function. We adopt the usual 
decomposition for the complex radius of curvature along 
its imaginary and complex part: 



1 



2i 
Kw 2 



Assuming that if 2(C) is real, and combining the imagi- 
nary part of Eq. (|A5|) with the real part of Eq. (|A6[) . one 
obtains: 



\A{Q\ 2 = \M 



irr. 



2 (0 



This relation reflects the conservation of the atomic flux 
T along the propagation. With our choice of normaliza- 
tion, the parameter |A| 2 is given by: 



(A7) 



3. Spherical Wave solutions to the nonlinear 
equation. 

With several approximations, it is possible to find sim- 
ilar solutions in the interacting case. Atomic interactions 
are described by the mean-field potential 

V l (x,y,C)=g I \^//(z)\ 2 H±(x,y,C)\ 2 with .9/° = ^^ 

which intervenes in the time-independent equation ver- 
ified by ip. Because we adopt here a different normali- 
sation for the wave-function, the nonlinear coupling con- 
stant g® differs from the coupling constant gj used previ- 
ously: <?J = 9i/N- The mean-field contribution induces 
the following transverse potential 

V±(x,y,Q = 5?I^//(C)| 2 (\^±(x,y,C)\ 2 ~ |^(0,0,C)| 2 ) 

in the paraxial equation verified by ip± . In the considered 
example, this potential receives no other contribution. 
The subsequent analysis requires three important ap- 
proximations. First, it uses the "aberrationless approx- 
imation" , which assumes that the wave-function follows 
the Gaussian profile (|A4[) in spite of the non-linearity. 
Second, it assumes that the transverse mean-field poten- 
tial is well-described by a second order expansion, 



V^x^O^^g^f/iO^-^-r 2 



(A8) 



The term G(() = ff/|^//(C)| 2 can be seen as the atom- 
optical equivalent of a third-order non-linear permittiv- 
ity. Third, it neglects the dependence on G(£) towards 
the altitude, which is a valid approach if the linear density 
n\D = mJ r /p(z) is a constant [xiii] . We assume from now 
on that the atomic flux T is constant and that the aver- 



and 



age longitudinal momentum p(z) = ■J 2m (E ~ V// (z)) 

Po// varies very slowly with z. The parameter £ can then 
be expressed simply as Q = m{z — z )/p //. Eq. 
the normalization of tp± [Eq. (|A7I) ] give readily: 



MO 



K irpo/fW^iC) 
Eq. (|A5j) can then be recast as: 
q'-l T ( 4 



T c \K 2 w\Q 







The quantity T Cl called critical flux, reads T c = 
irp //h 2 1 l (2g l jm 2 ). The last equation may be split into 
its real and imaginary part along: 



1 

R 2 



( 2 



V Kw 2 







(A9) 



[xiii] This approximation is indeed implicit in the treatment of 
Belanger and Pare |24| . since it is necessary to obtain the non- 
linear paraxial wave-equation which is the starting point of their 
analysis. 
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and 



/ 1 



V Kw 2 



R \Kw 2 







(A10) 



where we have introduced the dimensionless parameter 
a = 1 + T I T c . This system can be uncoupled thanks to 
the following trick: Eq. (|A10|) is multiplied by iy/a and 
added to Eq. (|A9[) . One obtains: 



/ 2 V 1 ( 2 



This equation can be simply interpreted as 
Qnl - 1 = 

with the generalized complex radius of curvature: 



qNL -{~ 



Kw 2 



Its very simple evolution 



Qnl (#) = qNL(za) 



m(z - go) 
PO// 



(All) 
(A12) 

(A13) 
(A14) 



gives readily the real radius of curvature R(z) and the 
width w(z) for any altitude z. One thus has, as in the lin- 
ear case, a simple spherical- wave solution (IA4|) . Indeed, 
this method allows one to approximate very efficiently 
the nonlinear propagation of a wave-function of initial 
Gaussian profile. Consider a Gaussian atomic beam of 
width w(zq) = wo at the waist (R(zq) = +oo) situated 
at the position z on the propagation axis. Eqs. (|A13[) 
and (|A14[) show that the beam width follows: 



w(z) 



r, 2 



a(z - z a ) 2 



(A15) 



The width of a self-interacting atomic beam evolves thus 
as an interaction-free beam in which the propagation 
length from the waist is multiplied by a factor ^fa. As 
far as the paraxial beam width evolution is concerned, 
self-interaction effects thus operate as a scaling transfor- 
mation of the free propagation with a factor i/a. The 
quantity i/a — 1 has the same sign as the scattering 
length a, so one checks that Eq. (|A15j) leads consistently 
to a faster expansion for repulsive interactions and to a 
slower expansion for attractive ones. As in optics, this 
treatment can thus be applied to discuss the self focusing 



for matter waves. It is, however, important to keep in 
mind its validity domain and the several hypothesis 
required - constant longitudinal velocity, cylindrical 
symmetry, paraxial propagation and Gaussian shape 
approximation -. Last , w e point out the independent 
work of Chen et. al. 58] on this nonlinear radius of 
curvature. 

APPENDIX B: SECOND-ORDER 
COMPUTATION OF THE NONLINEAR ABCD 
MATRIX. 

1. Expression of the second-order matrix. 

In this Appendix, we discuss the nonlinear cor- 
rections to the ABCD matrix associated with the 
second-order term of the Magnus expansion f2 2 (t,<o) = 
\ Jl dh dt 2 [N(ti), N(t 2 )], which reads: 



fi2(Mo) = 



f S(t,t ) 



-S(Mo) 



with S(t) = f dh f 1 dt 3 (7(ii) - 7(*a)) (Bl) 

This term, arising from the non-commutativity between 
the Hamiltonians taken at different times, naturally de- 
pends on the ordering chosen for the successive lenses. 
Because of the cloud expansion, lenses are ordered from 
the most divergent to the less divergent. To discuss 
the effect of this second-order contribution on the wave- 
function, it is useful to compute the exponential: 



exp 



[fi (2) (Mo)j 



,S(t,t ) 







e -s(t,t ) 



(B2) 



The action of such matrix onto the position-momentum 
width vector (X, Y), defined in Sec. IIII Ai would operate 
a squeezing between position and momentum. This 
squeezing is indeed a consequence of our aberrationless 
approximation, in which the propagation leaves the 
phase-space volume invariant: the expansion of the 
cloud size must be, in our treatment, compensated by a 
reduced momentum dispersion. One finds consistently 
that the diagonal matrix elements S xx ,yy,zz(t)-> involved 
in (|B2[) , are positive, which results from the decrease of 
the matrix elements "/ X x,yy,zz(t) with time. 

The ABCD matrix obtained from a second-order ap- 
proximation of the Magnus expansion reads: 



M M (tt x ,„n /cosh^(Mo) + ^(Mo)^ 7 ^T i C-'o) 2 ^ 
2 ( ' °' ° } - ^ Z { (7)^^ coshiMMo) " ^(Mo)^SS 



(B3) 



We have introduced the functions Ku(t,to) = yS^it^to) + (ju) (t — t ) 2 . An analytic expression of (7^) can 



13 



be found for cigar-shaped condensates in Eq. (|C1[) of Appendix [Cl The computation of the quantity Su(t,to) is 
straightforward, but it involves tedious algebra. Higher-order contributions to the ABCD matrix ([20]) involve various 
integrations which need to be performed numerically. 



2. Comparison with the first-order matrix. 



Let us expand the matrix (|B7|) in the short duration limit. We consider an atomic cloud initially described by a 
Gaussian wave- function (|10p of spherical symmetry i.e. w x0 = w y o — w z0 = wq. Such assumption does not change 
the nature of the discussion, but it considerably simplifies the algebra: the 3x3 matrices ^(t), S(t) and Kit) are 
then proportional to the matrix identity /3 and can be identified to scalars. j(t) can be expressed as a function of 
two time scales T\,t% involving the sample radius u>o, the scattering length a and fundamental constants: 



7(*) 



-5/2 



Ti 



mwQ 



-, 72 



Wq 

4:T:a 



(B4) 



The quantity S(t) (|B1[) can be expressed thanks to a second-order Taylor expansion of j(t). Setting r = t — t and 
noticing that 7' (to) — 0; one obtains: 



S(t) = 



5 r 4 



0(t 6 ) 



(B5) 



which yields for the quantity K{t): 



K(t) = V(l)r 1 + 



25 t 



72 T?T? 



O (r 8 ) 



(B6) 



Using this expansion and that of x 

( 



M^(t,X ) = M^'(t,X ) + 



(i), 



V 



sinhx/x, one can express the second-order matrix (r, Xq) as: 



T 4 sinh(( 7 ) 1 / 2 T) 



W 1 ^ 



cosh({ 7 ) 1 / 2 i 



dnh({ 7 ) 1 / 2 7 



5 T 4 S inh(( 7 ) 1 / 2 T) 
(7> 1 / 2 r 



0(r 8 0B7) 



This expansion shows that the first-order term is a valid approximation as long as: 



V47ra/ 



1/6 



(B8) 



Considering for an instance an initial cloud size of Wq = 25 /iin and the 87 Rb scattering length a = 5, 7nrn, one 
obtains t% — 0,14 s, ti = 1,63 s, and r c = 0,31s. Note that the relevant small parameter e, weighting the relative 
correction brought by the second-order term, decreases as e = (t/t c ) 4 when r/r c — > 0. 



APPENDIX C: ABCD MATRIX ELEMENTS FOR THE CIGAR-SHAPED CONDENSATE 



We evaluate in this appendix various primitives necessary to explicit the non-linear ABCD matrix to first 
order in the Magnus expansion given in Eq. (|22p . We consider a cigar-shaped cylindrical condensate with a long 
vertical extension: w x — w y — w r <^ w z . We remind the linear evolution of the width given by Eq. (1111) i.e. 

uv,zt = ^Jw^ z0 + Av^ z0 {t — t a ) 2 . We use the short-hand notation Aw rjZ0 = h/(mw r ^ z0 ). 
We seek to evaluate the average (7)^ = 1/r f t dtju(t) of the time-dependent coefficients: 
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with 70 = g//((27r) 3 / 2 m). These quantities are readily obtained: 



(7rr) = 7o 

(izz) = 70 



{Av 2 w 2 zQ - 2mg A^ ) Arctan A(t) ' 

2w 2 (Av%w 2 z0 - w 2 Q Av 2 zQ )w 2 f 2w%(Av 2 w 2 z0 - w 2 r0 Av 2 z0 ) 3 / 2 (t - t )_ 



A^ Arctan A(i) 



«4( A ^0 W rO - A «rO W ?o)^t ™ro(Av 2 Q W 2 zQ - W 2 Q Av 2 ) 3 / 2 (t - t ) 



with A(t) 



y/Av 2 Q w 2 z0 - w 2 rQ Av 2 z0 {t - t ) 



W rt W zt 
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